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Abstract. In this article we study the problem of thoracic image regis¬ 
tration, in particular the estimation of complex anatomical deformations 
associated with the breathing cycle. Using the intimate link between the 
Riemannian geometry of the space of diffeomorphisms and the space of 
densities, we develop an image registration framework that incorporates 
both the fundamental law of conservation of mass as well as spatially 
varying tissue compressibility properties. By exploiting the geometrical 
structure, the resulting algorithm is computationally efficient, yet widely 
general. 
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1 Introduction 

In this paper we consider the problem of tracking organs undergoing deformations 
as a result of breathing in the thorax and imaged via computed tomography 
(CT). This problem has wide scale medical applications, in particular radiation 
therapy of the lung where accurate estimation of organ deformations during 
treatment impacts dose calculation and treatment decisions [^[^[^[^. The 
current state-of-the-art radiation treatment planning involves the acquisition 
of a series of respiratory correlated CT (RCCT) images to build 4D (3 spatial 
and 1 temporal) treatment planning data sets. Fundamental to the processing 
and clinical use of these 4D data sets is the accurate estimation of registration 
maps that characterize the motion of organs at risk as well as the target tumor 
volumes. 

The 3D image produced from X-ray CT is an image of linear attenuation 
coefficients. The linear attenuation coefficient /i of a material is defined as 
/i = where am is the mass attenuation coefficient of the material and pm 

is the mass density. The linear attenuation coefficient is proportional to the true 
density and therefore exhibits conservation of mass. 

Currently, the application of diffeomorphisms in medical image registration is 
mostly limited to the L? image action of the diffeomorphism group, which is not 


a mass-preserving transformation. Furthermore, the diffeomorphisms estimated 
from typical image registrations algorithms (such as LDDMM or ANTS [^) 
do not accurately model the varying compressibility of different tissues types. In 
thoracic datasets, the lungs are highly compressible. Conversely, the bronchial 
tubes and the tissue surrounding the lungs are incompressible. During inhale, 
as air enters, the lung volume increases and the lung density decreases, while 
during exhale lung volume decreases and the lung density increases. But in both 
inhale and exhale, the lung mass is conserved. 

In this paper we use a cone-beam CT dataset of a rat acquired at 11 time 
points of an inhale-exhale breathing cycle. Figure shows the mass, volume, and 
density of the lungs of a rat at each time point of its breathing cycle, exemplifying 
these properties. 
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Fig. 1: Rat lung data: volume, density, and mass of the lungs during an inhale- 
exhale breathing cycle. As the volume increases, the density decreases, but mass 
is conserved. 


Both of these effects can be clearly seen in the histograms of a full-inhale and 
a full-exhale image, as shown in Figure 

In 2010, the EMPIREIO challenge compared registration algorithms 
applied to intra-patient thoracic CT images. The winner of the competition used 
an LDDMM method using normalized cross correlation metric 21 . This method 


does not model conservation of mass or spatially varying tissue compressibility. 
While others in this competition used the density action on these images [6|9| , none 
of these methods incorporate the spatially varying nature of tissue compressibility. 

We present an image registration technique that incorporates conservation 
of mass and organ compressibility. Instead of the image action of diffeomor¬ 
phisms, we use the physiologically appropriate density action. We also regularize 
the diffeomorphism by using a space-varying penalty which allows for high com¬ 
pressibility of the lung tissue while at the same time enforcing incompressibility 
of high density structures such as bone. The algorithm is based on the intimate 
link between the Riemannian geometry of the space of diffeomorphisms and the 
space of densities MPI . The resulting algorithm also has the added advantage 
that it is computationally efficient: orders of magnitude faster than existing 
diffeomorphic image registration algorithms. 











Fig. 2: Histograms of a full-inhale and full-exhale image. Each histogram has 
three peaks: the peak at 0 represents surrounding air, the middle peak represents 
lung tissue, and the peak at 90 represents soft tissue. For the lung tissue, the 
full inhale has higher volume but a lower image intensity than the full exhale, 
therefore showing conservation of mass. For the soft tissue, the average intensity 
does not change because it is incompressible. The slight drop in frequency of the 
full inhale is due to soft tissue leaving the image boundary. 


2 Mathematical Formulation 


Mathematically, the problem is to find a diffeomorphic (bijective and smooth) 
transformation between two densities on a subset i? C With a ‘density’ 
we mean a volume form on i?, i.e., an element of the form / dx where dx = 
dx^ A dx‘^ A dx^ is the standard volume element on and / = /(x) is a non¬ 
negative function on i?. The space of all densities on i? is denoted Dens(i?). One 
might, of course, identify I dx with its function /, and thereby think of Dens(i?) 
as the set of non-negative functions on i?. However, the invariance properties 
and geometry of the problem are remarkably more transparent when viewing 
Dens(M) as a space of volume forms. 

The group of diffeomorphisms Diff(i?) acts from the right on Dens(i?) by 
pullback: the action of ip e Diff(i?) on I dx e Dens(i?) is given by 

{(f, I dx) i-A (p*{I dx) = {\D(p\ I o ip) dx, (1) 


where \Dip\ denotes the Jacobian determinant of ip. The corresponding left action 
is given by pushforward: 

[ip, I dx) ip^{I dx) = ((/:?“^)*(/ dx) = (^\Dip~^\ I o ip~^) dx. (2) 


The Riemannian geometry of the group of diffeomorphisms endowed with a 
suitable Sobolev metric is intimately linked to the Riemannian geometry of 
the space densities with the Fisher-Rao metric. This has been developed and 
extensively studied in [4 13 14 : the basic observation is that there are Sobolev 














-metrics on the space of diffeomorphisms that descend to the Fisher-Rao 
metric on the space of densities. 

The distance associated with the Fisher-Rao metric is traditionally defined 
between probability densities (densities of total mass 1) and is given by 


= \/vol(i7) arccos 



( 3 ) 


where /io and jai are probability densities. It naturally extends to the space of all 
densities and the case when vol(i7) = oo, for which it is given by 


4(/o dx, h dx) = [ i^/To - y/hfdx . (4) 

JQ 

Notice that d?p{'^') in this case is the Hellinger distanee. For details, see [^. 

The Fisher-Rao metric is the unique Riemannian metric on the space of 
probability densities that is invariant under the action of the diffeomorphism 
group . This invariance property extends to the induced distance function, so 

d?p{lQ dx, Ii dx) = d‘^p{ip^{Io dx),(p^{Ii dx)) \/ip G Diff(i?) . (5) 



Fig. 3: Illustration of the geometry associated with the density matching prob¬ 
lem. The gradient flow on Diff(i?) descends to a gradient flow on the orbit 
Orb(/(ix, /o dx). While constrained to Orb(/dx, Jq dx) C Dens(i?) x Dens(i?), 
this flow strives to minimize the product Fisher-Rao distance to {{fo(p)dx, Ii dx). 


Motivated by the aforementioned properties, we develop a weighted diffeo- 
morphic matching algorithm for matching two density images. The algorithm 











is based on the Sobolev gradient flow on the space of diffeomorphisms that 
minimizes the energy functional 

E{ip) = dx), (/ o ip~^)dx) + dp{cp^{Io dx),Ii dx)). (6) 


This energy functional is only a slight modification of the energy functional studied 
in [^. Indeed, if / in the above equation is a constant a > 0, then (§1 reduces to 
the energy functional of Bauer, Joshi, and Modin §5.1]. Moreover, the geometry 
described in §5.3] is valid also for the functional (|^, and, consequently, the 
algorithm developed in § 5.2] can be used also for minimizing There the 
authors view the energy functional as a constrained minimization problem on 
the product space Dens(i7) x Dens(i7) equipped with the product distance, cf. 
Fig and § 5] for details on the resulting geometric picture. Related work 
on diffeomorphic density matching using the Fisher Rao metric can be found 



Using the invariance property of the Fisher-Rao metric and assuming infinite 
volume, the main optimization problem associated with the energy functional (§ 
is the following. 


Given densities /q dx, Ii dx, and f dx, find (p G Diff(i?) minimizing 


E{cp)= [ {^/\D^\-lffocp-^dx+ [ 
Jo Jo ^ 


Ii] dx . 


EiM 


E2i<p) 


( 7 ) 


The invariance of the Fisher-Rao distance can be seen with a simple change 
of variables x ^(^), dx \Dip\dy, and \D(p~^\ Then, Equation m 

becomes 

E{cp)= [ {1 - ^/\DP\ff dy + [ (^/^o-^y\Dcp\hop\y. (8) 

To better understand the energy functional E{}p) we consider the two terms 
separately. The first term E\[^') is a regularity measure for the transformation. 
It penalizes the deviation of the diffeomorphism (p from being volume preserving. 
The density f dx acts as a weighting on the domain i?. That is, change of volume 
(compression and expansion of the transformation (p) is penalized more in regions 
of i? where / is large. The second term E 2 {(p) penalizes dissimilarity between 
Jo dx and (p*{Iidx). It is the Fisher-Rao distance between the initial density 
Iq dx and the transformed target density dx). Because of the invariance (§ 

of the Fisher-Rao metric, this is the same as the Fisher-Rao distance between 
Ii dx and ip^{Io dx). 

Solutions to problem 0 are not unique. To see this, let Diff/(i7) denote the 
space of all diffeomorphisms preserving the volume form I dx: 


Diff/(i7) = {ipe Diff(i?) I \D(p\ {Io(p)= /}. 


(9) 









If (/:? is a minimizer of ^(•), then ^j; o ip ior any 

e Diffiyo(f?) := Diffi(f2) nDiff/o(f?) (10) 

is also a minimizer. Notice that this space is not trivial. For example, any 
diffeomorphism generated by a Nambu-Poisson vector field (see [^), with Iq as 
one of its Hamiltonians, will belong to it. A strategy to handle the degeneracy 
was developed in § 5]: the fact that the metric is descending with respect 
to the metric on Diff(i7) can be used to ensure that the gradient flow is 
infinitesimally optimal, i.e., always orthogonal to the null-space. We employ the 
same strategy in this paper. The corresponding geometric picture can be seen in 
Fig.[3l 


3 Gradient Flow Algorithm Development 


We now derive in detail the algorithm used to optimize the functional defined in 
Equation!^ The i7^-metric on the space of diffeomorphisms is defined using the 
Hodge laplacian on vector fields and is given by: 

Gl{U,V)= [ {-Au,v)dx . (11) 

JQ 

Due to its connections to information geometry we also refer to this metric as 
information metric. Let E denote the gradient with respect to the information 
metric defined above. Our approach to minimize the functional of (§ is to use a 
simple Euler integration of the discretization of the gradient flow: 


p = — E{p) 


( 12 ) 


The resulting final algorithm (Algorithm Q is order of magnitudes faster than 
LDDMM, since we are not required to time integrate the geodesic equations, as 


necessary in LDDMM 23 


In the following theorem we calculate the gradient of the energy functional: 
Theorem 1. The -gradient of the matching functional (|^ is given by 


= -z\-i( - v(/ o - Pd^\))- 

yiD</.-i|joo<^-iv(vT) + v(yiD</.-i|/oo<^-i 


(13) 


Remark 2. Notice that in the formula for E we never need to compute p, 
so in practice we only compute p~^. We update this directly via p~^{y) ^ 
p~^{y + E) for some step size e. 

Proof. We first calculate the variation of the energy functional. Therefore let ps 
be a family of diffeomorphisms parameterized by the real variable s, such that 


Ps =v op. 

s=0 


Pq = p and 


d 

ds 


(14) 







We use the following identity, as derived in 10 


ds 


=iy|D<y?|div(v) o (f. 


s=o ' ' ■ ' 2 

The variation of the first term of the energy functional is 
d 

= 

IQ 


(15) 


Ei{ip) = [ fix){^/\D(p{x)\ - l)^/\D(p{x)\div{v) o ip{x)dx (16) 
=0 JQ 

We do a change of variable x i—)■ dx i—?■ \Dip~^{y)\dy, using the fact that 

= \D<p-Hy)\'^ 


■ f ^(2/)(l - 'J\Dip i(2/)l)div(n)(2/)dy 

Jn 


{f°V> ^(1-V|£>¥’ 


= -(V 


1,2 (R3) 


(17) 

(18) 
(19) 


using the fact that the adjoint of the divergence is the negative gradient. For the 
second term of the energy functional, we expand the square 


E 2 {(p) = [ Io{x) - 2^Jlo{x)Ii O ip{x)\Dif{x)\ + ho (p(x)\D(p{x)\d 2 

Jn 


( 20 ) 


Now h o ip{x)\Dip{x)\dx is constant (conservation of mass), so we only need 
to minimize over the middle term. The derivative is then 


ds 


E2{(f) = - [ 2y^Io{x){V\^^v) O (p{x)y^\Dip{x)\ 

=0 Jn 

— \/7o(^)^^To~^(^)TD^(^div(n) o ip{x)dx. (21) 
We do the same change of variables as before: 

= -[ yjo ° , {2V,/T^^viy) + y7^div(n)(y)) (22) 

- 7o o (p-i/i,div(n)^^^^^^^ 

+ (v(yi7i^-Voo^-)vT,.)^,^^3^. 


(23) 


(24) 


From the above equations we conclude that: 

-A{V^‘e) = -V (/ o ^-1(1 - vlD^)) 

- v/|D^-i| Jo o (p-iV/fT + V (v/|Z)(p-i| 7o o ^-i) ^/h (25) 






























Since we are taking the Sobolev 
the right hand side of Equation 


^dient of E, we apply the inverse Laplacian to 


25 


to solve for E. 


Algorithm 1 Final Algorithm 

Choose e > 0 
Set = id 
Set \Dip~'^ \ — 1 
for iter — l...NumIters do 
Compute — E ^ 

Compute u = -V(/ o + V(V^*/o)VA 

Compute V — —A~^(u) 

Update 

Update \Dip~^ \ \Dip~'^ \ o 

end for 


Remark 3. Algorithm [I] constructs the mapping Lp~^ by numerically integrating 
the vector field v. Thus, for small enough e, the computed transformation (p 
a diffeomorphism (as is also the case in LDDMM). 


is 


4 Results 


We applied the proposed method to the previously mentioned rat dataset. In 
this dataset, an anesthetized rat was placed on a mechanical ventilator. This 
ventilator sent 11 gate signals to the cone-beam CT per breathing cycle, assuring 
that all projections would all be acquired at a consistent points of the breathing 
cycle 11 . Previous literature has shown that cone-beam CT is inadequate in 
estimating the true linear attenuating coefficient density [^, so we empirically 
estimated the density as the square of the the original data. 

For these results we estimated the deformation from the full-exhale to the 
full-inhale image. The deformation was computed on the resolution of the original 
3D volume (245 x 189 x 217); all the figures show the same 2D coronal slice of 
this volume. Shown in Fig. are the coronal sections of full exhale, the full exhale 
deformed via the density action, and the corresponding image at full inhale and 
the estimated deformation. 

For the compressibility penalty /, we used a soft thresholding of the intensity 
values of the initial image using the logistic function. High intensity regions of 
the CT image (corresponding to bone and soft tissue) were given a high penalty 
{f{x) = lOcr) and low intensity regions of the CT image (corresponding to air 
and lungs) were given a low penalty {f{x) = .la) (see Figure]^ 

We implemented the proposed algorithm and LDDMM on a single Titan-Z 
GPU (using the PyCA software package bitbucket.org/scicompanat/pyca ) 
for comparison. The difference images are pictured in Figure The problem of 
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Fig. 4: Density action results. This figure shows the lung image at the full exhale, 
the full exhale deformed via the density action, and the corresponding image at 
full inhale. Shown in the right panel is the estimated deformation. 



\Iin Iex\ \Iin *^*(4eccdx)| \Iin lex \ (LDDMM) 

Fig. 5: Absolute value of image differences: The left panel shows the difference 
between the original full exhale and the full inhale images. The center panel 
shows the result after registration using the proposed method. The right image 
shows the result using LDDMM with the image action. In LDDMM, there is 
significant error inside the lung due to the action not preserving mass. 


LDDMM using the action can be seen in this image. The Jacobian determinants 
are in Figure The proposed method constrains the contraction and expansion 
to inside the lung and outside the body. In this figure we also show the results of 
using the density action with a constant penalty function (/(x) = a.). 

The proposed algorithm is significantly faster than LDDMM; it runs at 400 
iterations per minute while LDDMM runs at 45 iterations per minute. We used 
10 time steps to integrate the geodesic equations associated with the LDDMM 
formulation. Since we are not required to integrate the geodesic equations in the 
proposed algorithm, we have nearly a lOx speedup compared to LDDMM. 


















Fig. 6: Jacobian determinants: On the left is the Jacobian determinant of the 
transformation estimated by the proposed method. Notice that the volume change 
is confined to inside the lungs and outside the body. In the center we use the 
density action, but without a local-varying penalty (i.e. f{x) = a). On the 
right is the Jacobian determinant using LDDMM. Without the local-varying 
penalty, there is contraction and expansion outside of the lungs. In LDDMM, the 
contraction and expansion outside of the lungs is even more severe. 




sig(x) 


Fig. 7: Energy plot and the logistic function used for the penalty. 


5 Discussion 

In this paper, we introduced a computationally efficient method for estimating 
registration maps between thoracic CT images. The proposed solution accurately 
incorporates the fundamental property of mass conservation and the spatially 
varying compressibility of thoracic anatomy. We conserve mass by viewing the 
images as densities and applying the density action of a diffeomorphism instead 
of the typical action. We limit the volume change in incompressible organs 
by placing a space-varying penalty on the Jacobian determinant of the diffeo¬ 
morphism. While any non-negative function f{x) can be used, we simply use 
a soft-thresholding function on the initial image. This choice is based on the 
assumption that low CT image values (such as the lungs and air) exhibit a large 




















amount of volume change whereas high images values (such as other soft tissue 
and bone) are quite incompressible. 
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